Spin waves and revised crystal structure of honeycomb iridate Na2lr03 
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We report inelastic neutron scattering measurements on Na2lr03, a candidate for the Kitaev spin 
model on the honeycomb lattice. We observe spin-wave excitations below 5 meV with a dispersion 
that can be accounted for by including substantial further-neighbor exchanges that stabilize zig-zag 
magnetic order. The onset of long-range magnetic order below Tn = 15.3 K is confirmed via the 
observation of oscillations in zero-field muon-spin rotation experiments. Combining single-crystal 
diffraction and density functional calculations we propose a revised crystal structure model with 
significant departures from the ideal 90° Ir-O-Ir bonds required for dominant Kitaev exchange. 

PACS numbers: 75.10.Jm, 78.70.Nx, 75.40.Gb, 61.72.Nn 



Transition metal oxides of the 5d group have recently 
attracted attention as candidates to exhibit novel elec- 
tronic ground states stabilized by the strong spin-orbit 
(SO) coupling, including topological band- or Mott- 
insulators [ij, quantum spin liquids [2|, field- induced 
topological order [s^ , topological superconductors [4] and 
spin-orbital Mott insulators j5|]. The compounds ^2lr03 
{A=Li, Na) d, 01) in which edge-sharing IrOg octahe- 
dra form a honeycomb lattice [see Fig. [Ua)], have been 
predicted to display novel magnetic states for composite 
spin-orbital moments coupled via frustrated exchanges. 
The exchange between neighboring Ir moments (called 



=1/2) is proposed to be [2j] 



JkS^S] + JiS.,-S. 



J' 



(1) 



where Jk > is an Ising ferromagnetic (FM) term aris- 
ing from superexchange via the Ir-O-Ir bond, and Ji > 
is the antiferromagnetic (AFM) Heisenberg exchange via 
direct Ir-Ir 5d overlap. Due to the strong spin-orbital 
admixture the Kitaev term Jk couples only the compo- 
nents in the direction 7, normal to the plane of the Ir-O-Ir 
bond [1,1^. Because of the orthogonal geometry, different 
spin components along the cubic axes (7 = x, y, z) of the 
IrOg octahedron are coupled for the three bonds emerg- 
ing out of each site in the honeycomb lattice. This leads 
to the strongly-frustrated Kitaev-Heisenberg (KH) model 
0, which has conventional Neel order [see Fig. [5^)] for 
large Ji, a stripy collinear AFM phase [see Fig. [3}:)] for 
0-4 ^ a ^ 0.8, where a = Jk/ (Jk + 2Ji) (exact ground 
state at a = 1/2), and a quantum spin liquid with Ma- 
jorana fermion excitations JLO] at large Jk (a > 0.8). In 
spite of many theoretical studies IllHlJI very few 
experimental results are available for ^2lr03 d, 0, [l5| . 
Evidence of unconventional magnetic order in Na2lr03 



came from resonant xray scattering jl5i] which showed 
magnetic Bragg peaks at wavevectors consistent with ei- 
ther an in-plane zig-zag or stripy order [see Figs. [5)d-c)]. 

Measurements of the spin excitations are very impor- 
tant to determine the overall energy scale and the rele- 
vant magnetic interactions, however because Ir is a strong 
neutron absorber inelastic neutron scattering (INS) ex- 
periments are very challenging. Using an optimized setup 
we here report the first observation of dispersive spin 
wave excitations of Ir moments via INS. We show that 
the dispersion can be quantitatively accounted for by in- 
cluding substantial further-neighbor in-plane exchanges, 
which in turn stabilize zig-zag order. To inform future ah 
initio studies of microscopic models of the interactions 
we combine single-crystal xray diffraction with density 
functional calculations to determine precisely the oxygen 
positions, which are key in mediating the exchange and 
controlling the spin-orbital admixture via crystal field ef- 
fects. We propose a revised crystal structure with much 
more symmetric IrOg octahedra, but with substantial de- 
partures from the ideal 90° Ir-O-Ir bonds required for 
dominant Kitaev exchange and with frequent struc- 
tural stacking faults. This differs from the currently- 
adopted model, used by several band-structure calcula- 
with asymmetrically-distorted IrOg octa- 
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tions 

hedra, with Ir-0 bonds differing in length by more than 
20%, improbably large in the absence of any Jahn- Teller 
interaction, and with the shortest Ir-0 bond length be- 
low 2 A, highly unlikely for a large ion like Ir''+. We show 
that the previously proposed structure is unstable with 
large unbalanced ionic forces, and when allowed to relax 
it converges to a higher-symmetry structure. 

As other "213" honeycomb oxides, Na2lr03 has an al- 
ternating stacking of hexagonal layers of edge-sharing 
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FIG. 1: (Color online) a) Layer stacking along the monoclinic 
c-axis with an in-plane offset along a (dashed box is the C2/m 
unit cell), b) Basal layer (z = 0) showing the Ir honeycomb 
lattice, c) Diagram to illustrate the layer stacking in the ideal 
honeycomb lattice. Ideal stacking of layers and stacking faults 
are explained in the text, d) Xray diffraction intensity in the 
(0, k, I) plane showing rods of diffuse scattering in between 
structural Bragg peaks along c* with selection rule h + k = 2n 
and k = 3m + 1 or 3m + 2 (n, m integers) modelled in e) 
by frequent in-plane translational stacking faults of the type 
shown by the thick arrows in c). 



NaOg octahedra and similar layers where two-thirds of 
Na are replaced by Ir to form a honeycomb lattice with 
Na in the center [see Fig.[T]D)]. To determine the precise 
structure xray diffraction was performed on a flux-grown 
single crystal of Na2lr03 0, Hal- The diffraction pattern 
showed sharp Bragg peaks which could be indexed by a 
monoclinic unit cell [see Fig. [T^)] derived from a parent 
rhombohedral structure with an ideal repeat every three 
layers. The monoclinic distortion leads to an in-plane 
shift of successive Ir honeycombs differing by 1.2% from 
the ideal value [— ccos/3 compared to a/3, see Fig. [1^)], 
well above our instrumental resolution, which enabled us 
to determine that our sample was a single monoclinic 
domain. The detailed refinement [l6| was performed us- 
ing both the published C2/c (No. 15) unit cell with 15 
refined atomic positions leading to values somewhat sim- 
ilar to Ref. [6] , and an alternative, higher-symmetry and 
half the unit cell volume, C2/m model (No. 12, shown in 
Figs. [T^-b) (as found for the related Li2lr03 [l^l), with 
only 7 refined atomic positions listed in Table HI Other 
structural motifs reported for "213" honeycomb oxides 
[lit including Na2Pt03, Li2Te03, Na2Tb03 were also 
tried but did not provide a good fit. We also tested for 
Ir/Na site admixture but this did not improve the agree- 
ment with data. 



TABLE I: Structural parameters extracted from single-crystal 
xray data at 300 K. (C2/m space group, a = 5.427(1) A, 
b = 9.395(1) A, c = 5.614(1) A, /3 = 109.037(18)°, Z=4). All 
sites are fully occupied. U is the isotropic displacement. The 
goodness-of-fit was 2.887 {Rtnt = 0.1247, = 0.0584) [l|j. 
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0.001(7) 



The C2/c structure can be described as a "super- 
cell" obtained from the C2/to structure by small dis- 
placements of atoms (of order a few % of the unit 
cell dimensions) leading to a doubled unit cell volume. 
Although C2/m and C2/c gave comparable agreement 
with the main Bragg peaks, the larger C2/c unit cell 
should be manifested experimentally by the appear- 
ance of new "superstructure" peaks at positions such as 
(odd, odd, half-integer) in the small unit cell description 
(C2/to). These superlattice peaks, however, were not ob- 
served in the data fli'], ruling out the C2/c model. Fur- 
thermore, in structural optimization calculations using 
VASPflp, [III (also confirmed by an all-electron LAPW 
code [20]) we find that the C2/c structural model, which 
has asymmetrically-distorted IrOg octahedra, is unsta- 
ble: (i) the forces on oxygen are very large, exceeding 3 
eV/A for the published C2/c cell |6j and (ii) when the 
structure is allowed to relax the oxygens move such as 
to recover the more symmetric C2/m structure with the 
Ir-0 distances converging to within 1.1% of the experi- 
mentally refined values in Table HI The IrOe octahedra 
are much more symmetric in the C2/m model with Ir-0 
distances and Ir-O-Ir bond angles ranging from 2.06 to 
2.08 A, and 98 to 99.4°, respectively, compared to the 
wider ranges 1.99 to 2.43 A , and 91 to 98° proposed 
before 0|. 

In addition to sharp Bragg peaks, visible diffuse "rods" 
of scattering were also observed [see Fig. [TJi)] and could 
be quantitatively understood [compare with calculation 
in Fig. [TJ;)] in terms of a structural model that allows for 
the possibility of faults in the stacking sequence along 
the c-axis. The stacking of atomic layers can be eas- 
ily visualized with reference to projections in the basal 
plane [Fig. [1];)], where A defines a nominal hexagonal 
lattice (made up of three triple-cell sublattices A1-A3), 
and B and C are also hexagonal lattices with positions in 
the center of a triangles of A sites. The atomic stacking 
is always in the ABC sequence to minimize the inter- 
layer Coulomb energy, i.e. Ir-O-Na-O-Ir-0 is Ai-B-C-A- 
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FIG. 2: (Color online) a) ZF /x+SR spectra on a polycrys- 
talline sample of Na2lr03 above and below Tjv. Solid lines 
are (top) guide to the eye and (bottom) a fit described in the 
text. b,c) Fitted parameters as a function of temperature. 



Bi-C. Only Ir layers have a sublattice index, indicating 
the position of the Na at the honeycomb center, as the 
other atomic layers are full hexagonal lattices. However, 
if neighboring Ir layers are only weakly interacting (as 
they are separated by a hexagonal Na02 layer) then the 
second Ir layer could be shifted to another position on 
the B-lattice, say B2 [thick arrows in Fig. [TJ;)] or B3, 
with only minimal energy cost, as that would not affect 
the bonding with the fully hexagonal Na02 layers below 
and above. To quantitatively verify this idea, we per- 
formed structural optimization calculations using VASP 
[1^ in an extended unit cell to include a stacking fault 
of the type illustrated in Fig. [!};) and found that the en- 
ergy cost of a stacking fault is extremely small, below 0.1 
meV/A^, explaining why such stacking faults are very 
likely to occur. 

The calculated scattering for such a microscopic model 
[lij indeed reproduces well the selection rule for where 
diffuse scattering occurs in Fig.[T}i-e). In particular there 
is no diffuse scattering along (OOZ), as this corresponds to 
adding all layers in phase irrespective of their in-plane 
translations. Also there is no diffuse scattering along 
(0,6n, {n integer), as again layers add in phase be- 
cause the two allowed in-plane translations have a phase 
factor equal to a multiple of 2tt. We use the strength of 
the diffuse scattering integrated between (020) and (021) 
relative to the intensity of the (020) peak (to have similar 
absorption factor), obtained experimentally as ~ 0.42, to 
estimate the probability for stacking faults p ~ 9%, this 
means that on average one fault occurs every 1/p ~ 10 
layers. We measured over 30 crystals from a batch and 
all showed diffuse scattering, suggesting that this is a 
common structural feature. 

Magnetic order of the Ir spins was detected by zero- 
field (ZF) muon-spin rotation (/i+SR) on a powder sam- 
ple of Na2lr03. Example raw spectra are shown in 
Fig. [2k)- At temperatures below Tjv = 15.3 K, we observe 
clear oscillations in the time-dependence of the muon po- 
larization, characteristic of quasi-static local magnetic 



fields at the muon stopping site. Fits to the time- 
dependent muon data reveal that two frequencies are 
present, indicating the presence of two distinct muon 
stopping sites with different local fields. The full spec- 
tra was fitted to the form A{t) — Aie^'^'i* cos(27rz/ii -|- 



<f)i) + ^26"-^^* cos(27rJ/2t + 02) + ^36"^* + Abg, where 
the last two terms account for muons polarized paral- 
lel to the local magnetic fields, and muons stopping in 
the sample holder (or cryostat tail), respectively. Using 
our best-fit parameters we estimate that the muons oc- 
cupy the two sites with a probability ratio of about 9:1. 
Both local fields set in at a common temperature, but 
have a distinctly different temperature dependence [see 
Fig. Eb)]. The relative weight of the second frequency 
component suggests that it may come from muon sites 
implanted near stacking fault planes, as such sites also 
occur in a similar proportion. Our value for Tn is consis- 
tent with both susceptibility measurements on the same 
batch, which indicated a clear anomaly (sharp downturn) 
near Tn as reported previously [3| , and the magnetic 
Bragg peaks observed in resonant xray scattering [l5| . 

The magnetic excitations were probed by powder in- 
elastic neutron scattering using the direct-geometry time- 
of-flight spectrometer MARI at ISIS with an optimised 
setup to minimise absorption [l^. Fig. [3^) shows the 
raw neutron scattering intensity as a function of wavevec- 
tor {Q — \Q\) and energy transfer deep in the ordered 
phase. An inelastic signal with a sinusoidal-like disper- 
sive boundary below a maximum near 5 meV is clearly 
observed at low Q. A gap, if present is smaller than 2 
meV. The magnetic character of the scattering is con- 
firmed by the broad, damped-out signal observed in the 
paramagnetic phase at 55 K [see Fig. [3t) and g) (con- 
trast filled and open symbols)]. Interestingly, the dis- 
persion boundary extrapolates at the lowest energies to 
a wavevector Q much smaller than that expected for 
conventional Neel order, Q(o20) = 1-34 A~^, so this 
magnetic order can be ruled out; in fact Q is close to 
the expected location of the first magnetic Bragg peak 
for both zig-zag or stripy order, (5(oio) — 0-67 A~^. 
Figs. [3]i) and i) show the calculated scattering from 
spin waves of a 2D Heisenberg model with up to 3rd 
neighbour exchanges, Ji,2,3, with zig-zag (Ji — 4.17 
meV, J2/J1 — 0.78, J3/J1 — 0.9) and stripy order 
(Ji = 10.89 meV, J2/J1 = 0.26, J3/J1 = -0.2), respec- 
tively (we neglect the interlayer couplings believed to be 
small). The constraints to reproduce the dispersion max- 
imum and the measured Curie- Weiss (CW) temperature 
(9 = -Sis + l)(Ji + 2J2 + J3)/kB ~ -125 K 0) are 
not sufficient to determine all 3 exchanges, so the values 
chosen are representative of the level of agreement that 
can be obtained [l^. The calculation for the zig-zag 
phase [Fig. [3h)] can reproduce well the observed disper- 
sion at low-Q (filled symbols), whereas the stripy phase 
[Fig. [3J)] cannot account for the strong low-Q dispersive 
signal and predicts stronger scattering at larger-Q's not 



4 



a) Neel 



b) zig-zag c) stripy 



d) O Neel 
^ zig-zag 
♦O stripy 




4 6 

Energy (meV) 



FIG. 3: (Color online) Diagram of a) Neel, b) zig-zag and c) 
stripy order, d) Reciprocal space diagram showing locations 
of magnetic Bragg peaks for various magnetic phases (inner 
hexagon shows first Brillouin zone of the honeycomb lattice), 
e) Powder inelastic neutron scattering data. The notable well- 
defined feature is the sharp lower boundary of the scattering 
at low Q (filled (magenta) symbols in h-j)), which we associate 
with a sinusoidal spin wave dispersion; this becomes damped 
out in the paramagnetic phase in f). Slanted thick dashed ar- 
row shows the scan direction in g). Gray shading marks the 
inaccessible region close to the elastic line dominated by inco- 
herent elastic scattering, g) Energy scan (solid points 4.6 K, 
open symbols 55 K) through the maximum spin-wave energy 
seen in e) fitted to a Gaussian peak (solid line), dashed line is 
estimated background, h-j) Calculated spherically-averaged 
spin-wave intensity [3| for the Ji,2,3 model with h) zig-zag 
or i) stripy order, and j) the KH model with stripy order for 
parameters given in the text. Solid red line in j) highlights 
the low-energy boundary, which coincides with the dispersion 
from r to the first softening point. 



seen. Calculations for the KH Hamiltonian ([T]) are shown 
in Fig. ISj) for a — 0.4 (lower limit for the stripy phase) 
and Ji = 25.85 meV to reproduce the CW temperature 
23 e = -S{S + l)(Ji - JK/3)/fcB. The lower bound- 



ary of the scattering at low Q (solid line) is predicted to 
have a quadratic shape near the first softening point, a ro- 
bust feature for any a throughout the stripy phase. This 
is in contrast to the data where the dispersion bound- 
ary (marked by filled symbols) has a distinctly different, 
sinusoidal-like shape with a curvature the opposite way. 
In addition, a different distribution of scattering weight 
to higher energies is predicted, but not seen in the data. 
We conclude that the KH model in the stripy phase has 
a qualitatively different spin- wave spectrum compared to 
the data. A minimal model that can reproduce the ob- 
served low-Q dispersion and which predicts distribution 
of magnetic scattering in broad overall agreement with 
the data up to some intensity modulations is shown in 
Fig. ^jp) and requires substantial couplings up to 3rd 
neighbors, which stabilize zig-zag magnetic order. Re- 
cent theory [isj proposed that in addition to couplings 
up to 3rd neighbors, a Kitaev term may also exist. We 
have compared the data with such a model as well [3] 
and estimate that a Kitaev term, if present, is smaller 
than an upper bound corresponding to a < 0.40(5). 

We note that sizeable Js's are not uncommon in trian- 
gular plane metal oxides. The reason is that even though 
Ji involves two hoppings and J3 four, the two additional 
hoppings are strong pdcr ones, and the hopping proceeds 
through intermediate unoccupied eg states [23] ■ In case of 
Na2lr03 the hopping proceeds through somewhat higher 
Na s orbitals, but these are very diffuse, and the corre- 
sponding tspcr parameter is sizeable. Near cancellation 
of the AFM and FM superexchange interaction for the 
nearest-neighbor path further reduces Ji compared to J3. 

To summarize, by combining single-crystal diffraction 
and LDA calculations we proposed a revised crystal 
structure for the spin-orbit coupled honeycomb antifer- 
romagnet Na2lr03 that highlights important departures 
from the ideal case where the Kitaev exchange domi- 
nates. We observed dispersive spin-wave excitations in 
inelastic neutron scattering and showed that substantial 
further-neighbor exchange couplings are required to ex- 
plain the observed dispersion and we proposed a model 
for the magnetic ground state that could support such a 
dispersion relation. 
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FIG. SI: (Color online) Radial distribution functions (RDFs) 
showing the difference in the local atomic environments for 
the previously reported C2/c structure [sj and the C2/m 
structure proposed in this study. In the C2/m case, the RDFs 
are plotted for the unit cell extracted from the experiment 
(red solid line), the unit cell with atomic positions relaxed in 
the DFT (blue dashed line), and the unit cell fully relaxed 
in the DFT (green dotted line). Note that a small Gaussian 
smearmg (cr = 0.008A) was used in the calculation of the 
RDFs. 



SI. Structural optimization calculations using 
VASP 



Supplemental Material for Spin waves and 
revised crystal structure of honeycomb iridate 
NaalrOa 

Here we provide additional information on 1) struc- 
tural optimization calculations to confirm the unit cell 
stability and estimate the energy of stacking faults, 2- 
3) the xray diffraction measurements and analysis of the 
diffuse scattering, 4) ^SR and 5) neutron scattering ex- 
periments, and 6-9) derive the spin-wave dispersion re- 
lations and dynamical structure factor in neutron scat- 
tering for the Heisenberg Ji.2,3, Kitaev-Heisenberg and 
Kitaev-Heisenberg- J2- J3 models for various magnetic or- 
ders. 



We used the Perdew-Burke-Ernzerhof (PBE) 
exchange-correlation functional Sl| within the gen- 
eralized gradient approximation (GGA) and the 
projector augmented waves method [S4I. The 2p semi- 
core electrons in Na were treated as valence. Numerical 
convergence was achieved with a 500 eV energy cutoff 
and dense Monkhorst-Pack fc-meshes [S3| of 7x7x3 for 
the previously reported p3 | C2 / c primitive unit cell and 
6x4x6 for the proposed C2/m conventional unit cell in 
Table I. We performed three types of calculations for 
the two structures: a static run with the experimental 
parameters, optimization of the atomic positions only, 
and full optimization of the atomic positions and lattice 
parameters. The residual forces and stresses were 
typically below 0.002 eV/A and 0.5 kbar, respectively. 
We found the magnetic and the spin-orbit interactions to 
have a rather small effect on the Na2lr03 structure and 
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the comparisons below are made for the non-magnetic 
case without the spin-orbit couphng. 

To iUustrate the differences in the local environments 
in Fig. lSll we plotted normalized radial distribution func- 
tions (RDFs) for all types of interatomic distances in the 
experimental and optimized structures. C2/c exhibits 
a considerable dispersion of the Ir-Ir and Ir-0 nearest 
neighbor distances critical for the magnetic ordering in 
the compound. The 0-Na and Na-Na separations are 
unphysically small and we observed large forces, over 6 
eV/A on Na and over 3 eV/A on O, at the beginning 
of the optimization run. The RDFs in C2/m with the 
experimental parameters demonstrate much more sym- 
metric local environments and a negligible variation of 
Ir-Ir lengths within the honeycomb lattice (below 0.3%). 
The calculated forces on atoms did not exceed 0.5 eV/A 
indicating a good agreement between the experiment and 
theory. Optimization of the atomic positions with fixed 
C2 / m experimental unit cell had little effect on the Ir-Ir 
distances because they are defined primarily by the in- 
plane lattice constants a and h. When fully optimized, 
C2/c and C2/m converged to the same structure with the 
C2/m space group and virtually indistinguishable RDFs. 
The enthalpy gains were 0.434 and 0.018 eV/atom, re- 
spectively (for comparison, the optimization of atomic 
positions in C2/m led to a 0.007 eV/atom gain). Note 
that the full optimization of C2 /m leads to ^ 2% elonga- 
tion of the Ir-Ir distances which is a typical bond overesti- 
mation observed for the GGA. For this reason we believe 
that use of the experimental lattice constants is more ap- 
propriate for the modelling of the magnetic interactions. 

To estimate the stacking fault energy we simulated 
Ixlxn [n — 2,...,6) supercells of the C2/m primitive 
12-atom unit cell with one Ir-Na layer and the two 
adjacent O layers shifted by 6/3 along [010]. The 
resulting lower-symmetry structures (C2 space group) 
had two stacking faults per unit cell and the same 
a X 6/2 — 25.49A^ x — y base. We optimized only 
the atomic positions keeping the experimental unit cell 
parameters fixed. The n = 2 structure gained additional 
symmetry operations (C2/c space group) upon relax- 
ation. The comparison of the faulted structures against 
the respective C2/m supercells with the same unit cell 
dimensions and the same A:-point meshes allowed us to 
reduce computational errors. However, the energy dif- 
ferences. En — EQ2/nn our uou- magnetic calculations 
without the spin-orbit coupling (SOC) proved to be 
exceptionally small in magnitude: 0.7, -1.7, -2.0, -2.6, 
-1.8 meV/(n x 12 atoms) for n — 2, ... ,6, respectively. 
For the smallest n ~ 2 structure we were able to 
calculate the energy difference with the FM ordering 
and the SOC as well and found En=2 — Ec2/m to remain 
small at 2.9 meV/(24 atoms). Based on these tests, we 
expect the stacking fault energy in C2/m to be below 
~ 0.1 me one to two orders of magnitude smaller 



k (r.l.u.) 
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FIG. S2: (Color online) Xray diflraction intensity in the 
(— 3, fc',/') plane: a) data, b) calculation for the "idealized" 
C2/c structure with atoms at special positions, equivalent to 
the C2/m model in Table I, c) calculation for the distorted 
C2/c model in [s3 | (assuming no Ir/Na site mixing in the hon- 
eycomb lr2/3Nai/3 layers). Notice the series of sharp peaks 
predicted in c) at (-3,1, even) positions, which however are 
not present in the data in a), d) Scan along the (-3,1,/') line 
(arrowed direction in a-c)) comparing data (solid circles) and 
calculation for the two models (triangles-C2/c and squares- 
C2/m). The calculated diffraction intensities have been mul- 
tiplied by an overall scale factor and have been convolved 
with a finite width Gaussian in momentum space to mimic 
the effects of the instrumental resolution. 



than typical stacking fault energies for elemental metals. 
For comparison, an ABCBA stacking fault generated 
by reflecting C2/m structure (which has the ABC ABC 
sequence along c) about a Na layer was calculated to 
have a much higher, measurable energy value of about 8 
meV/A^. 

S2. Xray diffraction and structural analysis 



X-ray diffraction was performed using a Mo-source 
Oxford Diffraction Supernova diffractometer on a single 
crystal of Na2lrOg of approximate size 220x 150x 10/im'^ 
grown via flux [Sj]. 96% out of over 1000 detected peaks 
were indexed by a single monoclinic domain. Structural 
refinement was performed using both a unit cell with 
space group C2/m, with parameters listed in Table I, 
as well as a unit cell with twice the volume and space 
group C2/c, using the SIR-92 and SHELX packages [S5 1 . 
The two unit cell parameters are related by a' = —a, 
b' = b, c' = a + 2c, c' = a? + ^c^ ^ 4accos/3, 
sin /?' = |f sin /?, and in terms of the reciprocal lattice 
components hi = — /i, fc' = —k, I' = h + 21, where primed 
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values refer to the C2/c model. Starting from the larger 
unit cell (C2/c) and slightly displacing the atoms to 
some "ideal" positions one recovers the higher-symmetry 
structure described by the smaller, C2/m, cell. The 
distinction between those two models is entirely due 
to such small atomic displacements, the presence of 
which is manifested in finite intensity diffraction peaks 
at {h',k',l') positions with h' odd, k' odd and I' even, 
which disappear when atoms are displaced to the "ideal" 
positions, when the structure recovers the C2/m sym- 
metry. This is illustrated by the calculated diffraction 
pattern in the (— 3,fc',/') plane where the "extra" peaks 
expected in the larger cell model C2/c shown in Fig. 
IS2b ) are not seen in the data plotted in Fig. \S2h ) , which 
is however fully consistent with the pattern expected for 
the higher-symmetry C2/m model shown in Fig. IS2b ). 
This is also apparent in Fig. IS2t i) showing a scan along 
the (—3, 1, /') line with extra peaks (triangles) predicted 
for I' = 0,2, not seen in the data (filled symbols). For 
completeness we note that we applied a shift of the 
fractional atomic coordinates in the C2/c unit cell (in 
the notation adopted in [SJ|) by (-1/4,-3/4,0) before 
converting them into fractional atomic coordinates of 
the C2/m cell (in the notation used in Table I), due 
to the different positions of the origin in the two space 
groups. 

53. Microscopic model of stacking faults 

The calculated diffraction pattern in Fig. le) was ob- 
tained numerically by direct structure-factor calculations 
using the DISCUS package |S6i]. We considered a "crys- 
tal" of 200a X 2006 x 4000c unit cells of NazIrOa (C2/m). 
To include the effect of stacking faults we assumed that 
each Ir layer has a choice with probability 1 — p to keep 
in-stacking-sequence with the layer below and p/2 to be 
shifted to either of the other two sublattice positions 
(translated in-plane by (0,1/3,0) or (1/2,1/6,0)), with 
p = for perfect stacking and p — 1/3 for a completely 
uncorrelated layer stacking sequence, a model first 
introduced to describe the stacking faults in the related 
material LiaMnOs [STj . 

54. Muon spin relELxation experiments 

Zero field (ZF) /i+SR measurements were made at the 
Swiss Muon Source (S/i+S), Paul Scherrer Institut, CH 
using the GPS spectrometer. For the measurement a 
250 mg powder sample of NaalrOs, which was used for 
inelastic neutron scattering measurement, was packed 
inside a silver foil packet (foil thickness 25 pm) and 
mounted on a silver sample holder. 

Fits of the data to an equation in main text reveal 
the evolution of Ui and Xi with temperature, as shown 
in Figs. 2(b-c). Unusually, the frequencies do not vary 
in fixed proportion, although they do tend to zero at the 



same temperature. The low-amplitude, higher frequency 
component V2 drops off far more dramatically than the 
large amplitude, lower frequency vi. In order to quantify 
this behavior, the frequencies were fitted to the phe- 
nomenological function i/,(T) = 1/^(0) [1 - (T/Tn)"']^'. 
A common value of Tn = 15.3(1) K was identified 
from fitting to this function. We find that a « 2 for 
both cases. The parameter /3 can be interpreted as an 
order parameter exponent. The other fit parameters are 
i/i(0) = 5.54(1) MHz, /3i = 0.36(1), ^2(0) = 6.20(3) MHz 
and /32 = 0.11(1). We note that Ai is an order of magni- 
tude larger than A2 , implying either that the distribution 
of fields is broader in the majority site or, assuming the 
fast fluctuation limit, that the fluctuation rate is smaller. 
The lower frequency oscillation, accounting for « 90% of 
the muon sites in the material, has a /3 value suggestive 
of the behavior of a three-dimensional (3D) system (for 
3D Heisenberg (3 = 0.367 and 3D Ising /3 = 0.326), while 
the minority muon site has an exponent value more 
similar to that expected for a 2D Ising system (for which 
/3 = 0.125). These seem to suggest that the magnetic 
fluctuations have a rather different character at the two 
muon sites. 

55. Inelastic neutron scattering experiments 

Inelastic neutron scattering measurements were made 
using the direct-geometry time-of-flight spectrometer 
MARI at ISIS using an incident neutron energy of 18 
meV, which covered the full bandwidth of magnetic 
excitations with a zone boundary energy near 5 meV. 
The instrumental energy resolution was 0.67(1) meV 
(FWHM) on the elastic fine. The sample was ^ 10 g of 
NaalrOa powder spread out in a very thin layer (< 1 mm 
to minimise neutron absorption) inside of an annular 
can of outer diameter of 40 mm and height 50 mm. 
Counting times for the data in Figs. 3e-f) were 28 and 
7 hours, respectively, at an average proton current of 
150/1 Amps. 

56. Spin-wave dispersions for the Heisenberg 
Ji_2,3 model in the zig-zag and stripy phases 

Here we outline the derivation of the linear spin wave 
dispersion relations and dynamical structure factors rele- 
vant for neutron scattering for various spin Hamiltonians 
on the honeycomb lattice. For the Heisenberg model with 
up to 3rd neighbour exchanges we extend previous results 
on the dispersion relations [S8| to include also the dynam- 
ical structure factors. For the Kitaev-Heisenberg model 
the spin- wave spectrum (including 1/S quantum correc- 
tions) has been studied before in a special "rotated" ref- 
erence frame [S9| . here we explicitly derive here the dis- 
persion relations and dynamical structure factors in the 
experimentally-relevant, un-rotated reference frame. For 
the Kitaev-Heisenberg- J2-J3 models both the dispersion 
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relations and dynamical structure factors have not been 
studied before. 

We start with the isotropic Heisenberg model on the 
honeycomb lattice with exchanges with up to 3rd nearest- 
neighbor, so called Ji,2,3 model with Hamiltonian 



INN 



2NN 



3NN 



where 1-, 2-, and 3NN indicate summing over all 1st, 2nd 
and 3rd nearest-neighbor pairs with couplings Ji , J2 and 
J3 [paths indicated in Fig. IS3b )]. where positive values 
correspond to antiferromagnetic exchanges. Depending 
on the relative ratio of the couplings there are six dis- 
tinct types of mean-field ground states Isi,|si^, which 
include the two candidate magnetic orders for Na2lr03, 
the zig-zag and stripy AFM orders shown in F igs. [S3k -b') 
(labelled II and IV, respectively, in iSSj, IsT^]). Both of 
those magnetic structures have four magnetic sublattices 
(labelled A-D) and can be described by a rectangular 
magnetic unit cell (dashed box in Figs. [S3^-b)), which 
coincides with the in-plane chemical unit cell a x 5 of 
Na2lr03. Within a single layer the Ir honeycomb lat- 
tice in very close to ideal {b/a ~ -s/S) in spite of the 3D 
monoclinic crystal structure, so we treat here the ideal 
2D honeycomb lattice with 3-fold symmetry. In this case 
the magnetic order can have three spacial domains, one 
such domain is shown for both structures in Figs. IS3b - 
b), the other two magnetic domains are obtained by ±60° 
rotation around the direction normal to the plane. 

Using a standard Holstein-Primakoff transformation in 
the large-S' limit the Hamiltonian becomes (to leading 
order) a quadratic form of magnon operators 



■H = ^ XtHX + 7V(1 + 1/S')£;_ 



MF 



(S2) 



where higher than quadratic terms are neglected. Here 
Em f is the mean- field ground state energy (per spin) and 
N is the total number of spin sites. The sum extends over 
all wavevectors q in the first magnetic Brillouin zone. 

For the zig-zag order in Fig. [S3h ) we define the opera- 
tor basis as X^' = f 



-q , b-q\ where a - 



- d label 



operators on sublattice A-D, i.e. (cq) creates (anni- 
hilates) a plane-wave magnon mode on sublattice A and 
so on. The Hamiltonian matrix in eq. (|S2[) is 



A B C D* 
B* A D C 
C D* A B 
D C B* A 



(S3) 



where 



A 
B 
C 
D 
V 



5{-Ji + 2J2- 
2SJir]cos{nh) 
2S'J2{cos [Tr{h - 

s {Jiv^ + J3 h 



3J3 + 2J2Cos(27r/i)} 



-k)] 
,-4 



-f cos [7r{h — k)]} 
27/2 cos(27r/i)]} 



Here {h, k) are components of the wavevector q in units 
of the reciprocal lattice of the ax b rectangular unit cell 
shown in Fig. [S3k ) . By periodicity the above expressions 
are valid for any momentum, not necessarily restricted to 
the 1st magnetic Brillouin zone. Di agon alisation of the 
Hamiltonian by standard techniques |Slll . S12] to obtain 
the normal magnon modes gives two doubly-degenerate 
dispersions 



A^ + BB* -C^- DD* 



±y/i\AB - CD*\^ - \B*D* - BD\ 



(S4) 



We have explicitly verified for the same model (|S1I) that 
eq. (|S4|) agrees with earlier results of [sl| [eq. (5.21)]. 
The spin-wave intensity in neutron scattering is propor- 
tional to the dynamical structure factor (expressed as 
S^^{Q,u}) for spin fiuctuations along the cc-direction and 
similarly for y-direction) and an analytical expression for 
this in the case of a Ham iltonian of the form in eq. (|S3[) 
are given explicitly |S12l | [eq. (A3)] and for brevity are 
not included here. 

The spin-wave dispersions in (IS4p (and their intensity 
dependence) for the zig-zag phase are plotted for repre- 
sentative exchange values in Fig. [S3h). As expected, the 
acoustic magnon, oj", is gapless with a linear dispersion 
at the magnetic Bragg peak at the Y point, is also linear 
and gapless at the X point, but has zero intensity because 
the structure factor for magnetic Bragg peaks also can- 
cels at this position. Furthermore, the dispersions soften 
and appear gapless at the M point and others part of 
the quartet (±1/2, ±1/2), which are Bragg peaks for the 
other two magnetic domains rotated by ±60°. This soft- 
ening is a general feature of linear spin-w ave d ispersions 
for a multi-domain magnetic ground state jS12| , however 
the fact that the gap closes at those points is not pro- 
tected by any symmetry, but is an artefact of the linear 
spin-wave approximation; by analogy with relat ed sp in- 
wave models for other multi-domain structures [S1^ we 
expect the dispersions to become gapped at the softening 
points when quantum fluctuations to 1st order in l/S* are 
included. 

A spherical average of the spin- wave spectrum (includ- 
ing various prefactors listed in eq. (|S9P below) is shown 
in Fig. 3h). The dominant contribution to the low-Q dis- 
persive edge of the strong signal near the first softening 
point ((5=0.67 A~^) is due to acoustic magnons on the 
Lo^ branch emerging out of the Y point and dispersing in 
the Y— )• r direction [see Fig. [S3li)] and also magnons on 
the branch emerging out of the M-point and dispers- 
ing in the M— >■ F direction. To reproduce the observed 
low-Q dispersion in the powder data we have imposed the 
constraint that the zone-boundary energy of the lowest 
branch on the F-Y line reproduces the observed maxi- 
mum of the low-Q dispersion, i.e. ui^ (O, 5) = 5 meV. 
This constraint together with the condition that the ex- 
changes reproduce the observed Curie- Weiss constant 
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a) zig-zag b) stripy c) reciprocal space x FY p' M F 




(-1 0) (0 0) (0 1) (1 1) (1/2 1/2) (0 0) (-1 0) (0 0) (0 1) (1 1) (1/2 1/2) (0 0) 



FIG. S3: (Color online) a) Zig-zag and b) stripy order. Dashed rectangular box of size a x b shows the magnetic unit cell 
containing four sublattices A-D. In a) solid, dashed and dot-dashed lines show paths for Ji, J2 and J3. b) Bond labels x, y, z refer 
to the components of the spins at the two bond ends coupled by the Kitaev term. Inset show projection of the (cubic) x, y and 
z axes onto the honeycomb plane, c) 2D reciprocal space showing magnetic Bragg peak positions for various magnetic orders, 
d-j) Spin-wave dispersions along symmetry directions in reciprocal space (arrowed path in c)) for the KH, Ji,2,3 and KII-J2-J3 
Hamiltonians for exchange values and magnetic orders listed in the legends. Wavevectors Q are expressed in reciprocal lattice 
units of the rectangular magnetic unit cell. Colour is the dynamical structure factor (convolved with a Gaussian in energy for 
visualization, full width at half maximum = 0.15Ji), isotropic for the Heisenberg model in g-h) {S^^{Q,ljj) = S^^ {Q,u))) and 
different for the two polarizations x,y for the KH model in d-f). i-j) Dynamical structure factor for the KII-J2-J3 model with 
the zig-zag structure in a), where the ordered moments are along a general direction x' in the xy plane and y' is a direction in 
this plane normal to x' . 



9 = ^S{S + l)(Ji + 2J2 + h)lkB = -125 K cannot 
determine all three exchange values Ji, J2 and J3, but 
allow for a one-dimensional family of solutions located 
on a curve in the parameter space {J2/ Ji, J'i/ Ji) (the 
dotted line in region II in Fig. [54]). All sets of exchange 
values part of this family are broadly consistent with the 
data. The level of agreement that can be obtained is il- 



lustrated in Fig. 3h) for one representative solution (red 
star in Fig. [S4| , chosen as it comes closest to reproducing 
also the intensity distribution at the lowest Q. 

We now turn to the alternative magnetic structure, the 
stripy order shown in Fig.[S3b). If the spin- wave operator 
basis is defined as = [a^ , 6^ , c_q , d-q] , then the 
Hamiltonian reduces to the same form as in eqs. (IS2IS3I) . 
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FIG. S4: (Color online) Classical phase diagram of the Heisen- 
berg Ji,2,3 model on the honeycomb lattice, eq. (|S1|I. showing 
the regions of stability for zig-zag (II) and stripy order (IV). 
Phase I is coUinear 2-sublattice Neel order, III and V are in- 
commensurate spiral phases, and solid lines are phase bound- 
aries .SIO ]. The dotted line inside region II indicates possible 
solutions for a minimal model to describe the spin dynamics 
in Na2lr03 obtained by imposing the constraints described in 
the text (the red star is a representative solution for which 
the full spectrum is shown in Fig. IS3h )V 



with magnon dispersions given by eq. (jS4p . but where the 
expressions for the A — D parameters are 



A 
B 
C 
D 
V 



{Ji + 2.h - 3 J3 - 
= S {Jur^ + J3 [ry* - 
= 2S'J2{cos [7r(/i-t-fc)] 
= 2SJir)-^ cosiirh) 



2J2Cos(27r/i)} 
277-2 cos(27r/i)]} 

-f cos [irih — k)]} 



The resulting spin-wave dispersions and intensities for 
representative exchange values are plotted in Fig. [S3k ). 
In contrast to the zig-zag phase, for the stripy phase the 
acoustic magnon, w~, is gapless, with a linear dispersion 
and finite intensity at both the X and Y points, as 
both are magnetic Bragg peaks with non-zero structure 
factor (X four times stronger intensity than Y). Again, 
due to the three domain structure there is softening 
of the dispersion with an artificial gapless point at M, 
which is expected to become gapped when quantum 
fluctuations beyond the linear spin-wave approximation 
are included, as discussed earlier. A spherical averaging 
of the spin-wave spectrum in Fig. IS3b ) is shown in Fig. 
3i), here the strongest signal at low energies is due 
to scattering from acoustic magnons near the X-point 
{Q = 1.16 A~^) with weaker scattering from magnons 
near Y {Q — 0.67 A^^) and intensity decreasing rapidly 
for magnons with smaller momentum. 



S7. Spin-vi^ave dispersions for the Kitaev- 
Heisenberg model in the stripy phase 

For the nearest-neighbor Kitaev-Heisenberg (KH) 
model in eg. (1) the stripy phase in Fig. IS3b ) is the 
stable ground state for 0.4 < a < 0.8, where a = 
Jk/ {Jk + 2Ji) jS9|. This ground state is exact at a = 
0.5, when upon rotation of the coordinate system at 
certain sites the Hamiltonian (1) converts to that of a 
Heisenberg ferromagnet in a rotated basis [so'] . 

For each of the three bonds coming out of a honey- 
comb lattice site the Kitaev term Jk couples different 
spin components x, y, z expressed in terms of an orthog- 
onal (cubic) reference frame. This is oriented with the 
cubic [111] axis normal to the honeycomb plane and the 
projections of the i, y and z axes in the plane making 
120° as shown in Fig. IS3b ) inset. Each bond is labelled 
with the type of the spin component for the moments at 
the two bond ends coupled by the Kitaev term, i.e. the 
z-bond AB stands for exchange —JkS\S^ and x-bond 
AD stands for — JxS'^S'f, and so on. 

Due to the anisotropic nature of the Kitaev ex- 
change more coupling terms between magnon opera- 
tors on the 4 different magnetic A - D sublattices are 
generated as compared to the Heisenberg Ji,2,3 model. 
Thus, one needs to use the full 8-term operator ba- 



sis 



-q ' 



-q ' 



<3 ' 



, for 



which the Hamiltonian expressed in magnon operators to 
leading order still has the quadratic form (IS2p with the 
matrix H given by 
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(S5) 



where 



A =S{.h + JK) 

B =SJiri-^ 

C = —SJy^i sin (tt/i) rj 

D = S{2Ji - Jk) cos (tt/i) 7? 



Diagonalization to get the normal magnon modes jSll | 
gives four dispersion relations 

ujl^{q) = -DD* + \B-C\'^ 

±^AA^\B - C|2 - \D*{B - C) - D{B* - C*)|2 

ul^{q)^A^-DD* + \C + B\' 



±y/4A^\B + C|2 - \D*{B + C) - D{B* + C*)|^ 



(S6) 
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The dispersion curves are plotted for a = 0.4 in Fig. 
[S3H ) and a = 0.5 in Figs. lSSl e-f). where the colour repre- 
sents the dynamical structure factor, plotted separately 
for the spin fluctuations along x and y-axes, the presence 
of Kitaev bond directional exchanges make those the dy- 
namical structure factor non-equivalent. The structure 
factors were obtained from the eigenvectors of the Hamil- 
tonian matrix H in eq. (jSSP , using a numerical implemen- 
tation of a general algorithm to diagonalize a quadratic 
form of boson operators proposed in S14l | . Changing the 
relative strength of the Kitaev term, for example a = 0.4 
compared to 0.5, does not change the spectrum qualita- 
tively only introduces a weak dispersion in the gapped 
uji^2 modes, compare Figs. [S3h-e). 



The dispersions show many distinct features compared 
to the case when the same stripy ground state was stabi- 
lized instead by isotropic Heisenberg exchanges shown in 
Fig. IS3^). Notably there is no longer a gapless mode at 
the r point and at the Bragg peak positions (X and Y). 
The lowest mode softens at the M point as in previous 
cases due to the 3-domain structure of the stripy ground 
state. The dispersion is gapless at this point in the linear 
spin- wave approximation and a gap is predicted to open 
up when quantum fluctuations to 1st order in 1/S are 
included for any general a, except for the exactly solv- 
able point a = 0.5 where due to an exact cancellation 
the spectrum is gapless [S9| . 

A spherical average of the spin-wave spectrum in Fig. 
[S3H ) (including both the and S^^ dynamical struc- 
ture factors) is shown in Fig. 3j), the lower boundary 
of the scattering at low-Q (emphasized by the red solid 
line) is due to scattering off magnons on the 104 F-M 
dispersion branch near the M point. 

S8. Spin-wave dispersions for the Kitaev- 
Heisenberg- J2- J3 model in the zig-zag phase 



Here we explore the effects of adding a small Kitaev 
interaction Jk to the Ji,2,3 Hamiltonian when the ground 
state order is the zig-zag phase (this has rec ently been 
shown to be stable for a range of Jk values [S15| ). We 
obtain the spin- wave Hamiltonian matrix in this case by 
combing eqs. ((S3|) and ((S5)) as 
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(S7) 



where 

A 
B 
C 
D 
E 
F 

c 



S {- Ji + 2 J2 + 3 J3 + 2 J2 cos(27r/i) + Jk} 

{-1/2)SJkV-^ 

S {2JiCOs(7r/i) - (1/2)JkC}?7 

S { Jir;-2 + J3 [7,4 ^ 27;-2 cos(27r/i)] - Jk^?" 72} 

2S'J2{cos [7r(/i + k)] + cos [nih - k)]} 

{1/2)SJk(:v 

gfc7ri/3 



Diagonalization leads to four dispersion branches 



W?2 



(q) ^A^ 



E^ 



\D-F[ 



^3,4 



±y^i\A{B-C)+E{D- 
(q) =A^-E^ + \B + C\''- 



\D + F\^ 

± V4|A(B + C) - E{D + F)\^ - 52 



(S8) 



where 6 ^ \{B - C){D* - F*) - [B* - C*){D - F)\. 
The dispersions are plotted in Figs. IS31-i) for a^O.4 
(Jk/Ji = 4/3), J2/J1 =0.23 and J3/J1 =0.51. To 
discuss the key features of the spectrum it is helpful to 
visualize the degeneracies associated with the magnetic 
order. The magnetic structure is the zig-zag pattern 
shown in Fig. IS3h ) but where the spin direction can be 
either along the x direction to satisfy the Kitaev term on 
the x-type AD bond, or along the y direction to satisfy 
the Kitaev exchange on the y-type BC bond. At the 
classical level any in-between direction, i.e. in the xy 
plane, also has the same energy, so one expects a gapless 
mode associated with rotations in this "easy" plane. 
Indeed Fig. IS3j ) shows that the dispersion is gapless at 
the Y point and with strong intensity for fluctuations 
in this easy-plane (along the y' normal to the ordered 
direction labelled x'), and gapped for fluctuations along 
z out of the easy plane, see Fig. IS3i ). Furthermore, 
due to the honeycomb lattice geometry the magnetic 
structure is degenerate with another two domains 
rotated by ±60° around the axis normal to the plane, 
so the spectrum is gapless at the Bragg peak positions 
of those other two domains, at points equivalent to M. 
The Hamiltonian however does not poses any continuous 
rotational symmetry in the presence of the Kitaev term, 
so one might expect that small gaps would open at 
both Y and M points when quantum fluctuations are 
included so the spectrum would be fully gapped. For 
completeness we quote the Curie- Weiss temperature 
for this model 9cw = -S{S+l){Ji+2J2 + J3-JK/i)/kB- 

S9. Spherically-averaged neutron scattering 
intensity 

The one-magnon neutron scattering intensity including 
the magnetic form factor and neutron polarization factor 
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is proportional to 



(f/(g) 



1 - 



(S9) 

where we u sed for f{Q) the Ir*+ spherical magnetic form 
factor [S16j and assumed the factor equal to 2. Here 
Qx iQy) ^-re the components of the wavevector transfer Q 
along the x-axis (y-axis), where z is the ordered spin di- 
rection. The precise direction of the ordered moments 
(£-axis) with respect to the crystallographic axes has 
only a small effect on the powder-averaged spectrum via 
small intensity modulations through the polarization fac- 
tors ^1 — ^5^^ , however for concreteness, we included a 
specific moment direction for the comparison with data. 
For the Ji,2,3 model in Figs. 3(h-i) the moments were as- 
sumed to be aligned along the crystallographic a-axis (as 
suggested by resonant xray data [S17| ) and for the KH 
model [Fig. 3j)] the moment is assumed to be along the 
cubic f -axis closest to the a-axis (tilted out-of-plane by 
35.26° from the -a axis, see Fig. [S3b) inset). Eq. ((S9)) 
was numerically averaged over a spherical distribution of 
orientations for the wavevector transfer Q and convolved 
with the instrumental resolution to obtain the plots in 
Figs. 3h-j), directly comparable with the raw neutron 
scattering data in Fig. 3e). For the KH-J2-J3 model the 
intensity is also given by eq. (|S9[) but with the axis labels 
{x,y,z) replaced by (?/', z, x'), where the a;'-axis defines 
the ordered moment direction (located in the original xy 
plane) and y' and z are orthogonal directions to it. 
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